Method for Protein Structure Alignment Using Cellular Automata

ABSTRACT

The detection of protein sequence alignments firstly generates protein blocks having n successive C-alpha atoms from two protein sequences. A matrix of the distance of protein atoms from said C-alpha atoms for each block is then determined. A difference matrix is determined from respective distance matrices representing the differences between the protein blocks of the two protein structures. A Cellular Automaton (CA) model is generated in an initial configuration based upon the difference matrix, and evolved through at least one generation using predetermined rules. The rules determine which CA cells remain live and which are set to be dead.

CROSS-REFERENCE TO RELATED APPLICATION

This application is a divisional of prior application Ser. No. 10/940,192, which was originally filed on Sep. 14, 2004.

FIELD OF THE INVENTION

The present invention relates to protein structure alignment using cellular automata.

BACKGROUND

At the end of the year 2003 more than 23,000 protein structures had been deposited in the Protein Data Bank (PDB), with more structures being discovered with each passing day. This surge in structure information has lead researchers to look out for efficient techniques to compare protein structures, detect motifs, classify proteins under specific families, etc.

While the current dogma of genetics connecting sequences to structures (sequence→structure→function) suggests that it would suffice if sequences were studied in greater details to detect similarity in proteins and classify them, there are several instances when different sequences yield the same structure. Hence there is a concerted effort to work with the 3-dimensional structure of proteins directly.

With all this structural information overflow emerges new necessities: that of identifying similar structures and mapping them to families, a quick and fast way to detect similarity, identify motifs, find longest contiguous alignments, etc.

Distance matrices are known to be used for various protein structure-related work. DALI (proposed by Holm and Sander) is a well-known structure alignment algorithm utilizing the concept of distance matrices. In DALI, the three-dimensional coordinates of each protein are used to calculate residue-residue (C_(alpha)-C_(alpha)) distance matrices. The distance matrices are first decomposed into elementary contact patterns, e.g., hexapeptide-hexapeptide submatrices. Then, similar contact patterns in the two matrices are paired and combined into larger consistent sets of pairs. A Monte Carlo procedure is used to optimise a similarity score defined in terms of equivalent intramolecular distances. Several alignments are optimised in parallel, leading to simultaneous detection of the best, second-best and so on solutions.

A need exists, however, for an improved manner of processing protein structure information.

SUMMARY

A protein alignment detection system uses an algorithm incorporating Cellular Automaton (CA) models for obtaining pairwise structure alignment of two proteins. Two-dimensional CA models are used, and the alignment information is presented in the form of patterns formed by the cells of the CA.

Breaking down the protein structure into distance matrices of 5 peptide units (Penta Peptide Blocks), the described system uses the differences of these matrices to construct the 2-dimensional CA grid. Each cell in the CA grid corresponds to an alignment with the state of the cell (dead or alive) indicating the presence or absence of an alignment at that location as determined from the state of the surrounding cells. Starting from an initial unaligned state, the CA evolves through several generations according to a defined set of local rules. The emergent pattern, as the CA evolves through successive generations, yields the alignment.

The described system exhibits “emergent” behaviour. Each cell behaves in a strictly microscopic way, but each individual cell's behaviour leads to a macroscopic long range behaviour exhibited by the entire system, which collectively indicates the alignment.

DESCRIPTION OF DRAWINGS

FIG. 1 is a schematic representation of a protein alignment detection system incorporating cellular automaton (CA) models.

FIG. 2 is a schematic of a protein structure.

FIGS. 3A and 3B are the distance matrices for two protein structures.

FIG. 4 is the difference matrix computed from the PPB of the two proteins of FIG. 3A and 3B.

FIG. 5 shows an elementary CA iterated over one time step.

FIGS. 6A-6D show the components of a 2-dimensional CA.

FIGS. 7A-7D show a schematic representation of the generation of a cellular automation (CA) from a difference matrix.

FIGS. 8A and 8B show two sample protein structures.

FIG. 9 shows a difference matrix for the protein structures of FIGS. 8A and 8B.

FIG. 10A and 10B is a schematic representation of neighboring units for which the centre unit lives, and FIG. 10B is a representation of neighboring unit configurations for which the central unit dies.

FIG. 11 shows an initial protein structure configuration and 4-generation iteration of a difference matrix.

FIGS. 12-15 show the iterative protein structure over 3 generations for different protein-pairs.

FIG. 16 is a schematic representation of a computer system suitable for performing the techniques described herein.

DETAILED DESCRIPTION Overview

FIG. 1 schematically represents the operation of the protein alignment detection system described herein. Two protein structures P1 10 and P2 10′, which are to be aligned, are both deconstructed to smaller units [P1₁, P1₂, . . . and P2₁, P2₂, . . . ] by a respective PPBs generator 20, 20′. Each of the protein sub-units (P1, P2) comprise n successive C-alpha atoms of the primary amino acid sequence. An appropriate value of n is selected. The PPB generators 20, 20′ then compute the distance matrix of each PPB by determining the distance between protein atoms from the atom's C-alpha position.

For the given two protein structures, a Difference Matrix (DMX) is computed from the respective distance matrices by a generator 30.

A CA generator 40 uses the DMX to generate an initial configuration of a 2-dimensional CA. The CA model evolves through several generations according to a pre-defined set of microscopic rules. The rules are designed so that only cells corresponding to the longest alignments of the two structures survive.

Selection of Number of C-atoms (n)

Distance Matrices are a suitable technique for structure comparisons and alignments. There are a range of values that can be used for the protein sub-units. Choosing a value of “1” would not serve the purpose, because that would mean that the algorithm is dealing with individual points and there is no concept of a distance, bringing with it the problems of rotation, translation, etc. Any n>1 would be appropriate. On the other hand, choosing too high a value would mean that if there are motifs comprising of, say, 10 amino acid sequences, they would get lost in the resultant huge matrix since the difference matrix generation wouldn't give a low difference. Small structural details like loops and bends can be detected if n<10. The value of n is selected as being 5.

Distance Matrix

The Penta Protein Block (PPB) of a protein structure is denoted by Pi_(j), where i is the protein index and j is the PPB number. The set of PPBs for a protein structure P1 shown in FIG. 2 is:

-   -   PPB1=[A1, A2, A3, A4, A5]     -   PPB2=[A3, A4, A5, A6, A7]     -   PPB3=. . .

The protein structure P2 is not shown.

If A1 and A2 are the amino acids of a protein structure with their C-alpha positions given as A1 (X₁, Y₁, Z₁) and A2 (X₂, Y₂, Z₂), then the distance between them is constructed as:

D ₁₂=√((x ₂ −x ₁)²+(y ₂ −y ₁)²+(Z ₂ −Z ₁)²)   [1]

The distance matrix P1₁ for protein structure P1 is shown in FIG. 3A. The distance matrix P2₁ for protein structure P2 is shown in FIG. 3B. Both have been calculated by equation [1].

Difference Matrix

The difference matrix DMX is a (m,n) matrix given by:

DMX _(ij) =Σ|P1_(i) −P2_(j)| for i=1 . . . m and j=1 . . . n   [2]

FIG. 4 provides an explanation of equation [2] with an example of how the first cell of the DMX is calculated. DMX₁₁ is the sum of the absolute differences of each element of P1₁ and P2₁:

DMX ₁₁=|0.0-0.0|+|1.3-1.1|+|2.4-2.2|+|3.2-4.2|+|4.2-5.2| . . . |4.2-|5.2|+|3.6-4.6|+|2.1-3.1|+|0.9-1.2|+|0.0-0.0|=15.4

Similarly, DMX₁₂ would be calculated using P1₁ and P2₂, and so on, till all elements of the DMX are calculated. The completed DMX provides the basis for generating the initial configuration of the Cellular Automata.

Cellular Automata

Cellular Automata (CA) are a simple class of systems that exhibit the following characteristics.

-   -   Large number of homogenous components (simple finite state         machines)     -   Extended in space     -   Absence of central control     -   Limited communication among components.

A one-dimensional CA consists of a lattice of N identical finite state machines (cells) each with an identical topology of local connections to other cells for input and output. Each cell is indexed by its site number i=0,1,2 . . . N-1. A cell's state at time t is give by s^(t) _(i). The state s^(t) _(i) of cell i together with the states of the cells to which it is connected is called the neighborhood n^(t) _(i) of the cell i. Each cell obeys the same transition rule Φ, that gives the updated state s^(t+1) _(i)=Φ(n^(t) _(i)) for cell I as a function of n^(t) _(i). For a 1-dimensional CA, the rule table is:

TABLE 1 Neighborhood 000 001 010 011 100 101 110 111 η: Output bit Φ 0 1 1 1 0 1 1 0 (η):

FIG. 5 shows the components of a one-dimensional, binary state, r=1, elementary CA shown iterated one time step.

A 2-dimensional CA has a grid of cells as opposed to a single row of cells as in the case of a 1-dimensional CA. FIGS. 6A-6D shows the components of a 2-dimensional CA in greater detail.

FIG. 6A shows the initial configuration (IC) of the CA. The IC is the entire population of cells (both dead and alive) with which the CA starts evolving according to the defined rules. Cells may be dead (indicated by the value 0) or alive (indicated by the value 1). Cells may have other values in case of non-binary CA. In binary CA, however, “0” and “1” are the only permitted values.

FIG. 6B shows one of the rules of the CA. This rule says that if the c1, c3, c5 and the c7 neighbors of a cell are alive in any generation, then this cell lives in the next generation. There will be many such rules indicating the various possibilities of the neighbors' states and the states of the cell itself. This is an illustrative example showing only one such rule. In fact, there will be a total of 2⁸=256 such rules.

As shown in FIG. 6C, each cell has 8 neighbors marked 1-8 (this is also called the 9-neighbor square). The cell itself is marked 0. In light of this, the previous rule states that if cells 1, 3, 5 and 7 are alive, then the cell marked 0 is alive in the next generation.

Finally in FIG. 6D, the CA evolves from its initial configuration to its first generation. Shown is a cell which is dead in the initial configuration. Now, when the CA evolves into the next generation, the update rules are checked. One such rule says that if the 1, 3, 5, 7 cells are alive, then that cell lives. In this example, this rule is satisfied and so the cell becomes alive in the subsequent generation, as shown.

Generation of Cellular Automata

The Difference Matrix generator 30 takes the difference matrix and using the values of the matrix, transforms it into the initial configuration (IC) of the CA.

FIG. 7A shows the sample DMX. Each element in the DMX acts as a cell. The neighbors of each cell are examined to find the minimum among them and these neighbors are made alive. Consider the cell with value 3.8. The first grid in FIG. 7B shows the neighbors of the ‘3.8 cell’ and it can be seen that the minimum value neighbor is that with a value of 0.2. The second grid in FIG. 7B shows the live cell and its place in the overall IC.

This process is performed for each cell within the DMX to give the complete IC. The grid in FIG. 7C shows the completed IC, whereas FIG. 7D shows another representation, where dead cells (value-0) are marked black and live cells (value-1) marked white.

Alignment Principle

All alignments (significant or otherwise) between the two protein structures will be visible in the Difference Matrix. Patterns thus need to be identified.

FIGS. 8A and 8B are protein structures, and the shaded portions of the protein structures are evidently similar. The PPBs of protein structure 1 are:

P1₁ =A1-A5, P1₂ =A3-A7, P1₃ =A5-A9, P1₄ =a7-A11, P1₅ =A9-A13, P1₆ =A11-A15, P1₇ =A13-A17, P1₈ =A17-A21, . . .

The PPBs of protein structure 2 are:

P2₁ =A1-A5, P2₂ =A3-A7, P2₃ =A5-A9, P2₄ =A7-A11, P2₅-A9-A13,P2₆ =A11-A15, P2₇ =A13-A17, P2₈ =A15-A19, . . .

The DMX is computed and is shown in FIG. 9. The similar structures correspond to:

-   -   A7 to A18 of protein structure 1 (corresponding to P1₄ to P 1₇         PPBs), and     -   A3 to A15 of protein structure 2 (corresponding to P2₂ to P2₅         PPBs).

Looking at the DMX, it is evident that there will be a contiguous run of diagonal cells of low value corresponding to the alignment because the distance matrices will be similar at these locations. The location of these low valued cells will be corresponding to the locations of P1₄ to P1₇ against P2₂ to P2₅.

Notice from the DMX in FIG. 9 that alignment patterns can be formed as cross-diagonals running in either direction (i.e. either in “3-7” or the “1-5” directions, as defined in FIG. 6C).

It is evident in FIG. 9 that the long diagonal running in 3-7 direction comprises an alignment. The objectives of alignment are to:

-   -   Eliminate all spurious alignments;     -   Retain the largest contiguous alignments; and     -   Indicate the direction of the alignments.

CA Rules for Continuous Alignment

Detecting the significant contiguous alignments involves designing a set of rules using which, in use, causes the CA to evolve and in the process eliminate all other cells that are not part of the alignment. Instead of looking at all the eight neighbors of a cell for writing the rules, only the cross neighbors are of importance. The guiding rule is that a cell that is alive should remain alive into the next generation if and only if cross neighbors along one or both directions are alive, i.e. either the “1-5” or the “3-7” direction (refer again to FIG. 6C for the definition of these directions). Applying these rules successively to all the cells of a particular generation should eliminate all such cells that are not part of the alignment. FIG. 10A shows all possible cross-neighbor combinations for which the centre cell lives, while FIG. 10B shows the remaining cells with states not found in FIG. 10A that will die.

FIG. 11 shows a CA initial configuration of the “1a4f_A” and “1c40_A” protein structures and the same after the CA has evolved for 4 generations ((a)-(d)). As the CA evolved through several generations starting from the initial configuration, only those cells that contribute to an alignment survive and the rest of the cells are dead.

It can be observed that only those cells contributing to the alignment of similar structures survive in the final generations. At the same time, the live cells at the ends of the largest alignment are successively getting eliminated. However, this is not a cause of concern because the survival of a lone live cell (a cell without any live neighbors) in the nth generation implies that it is part of an alignment that is at least of length 2n+1 (n cells on either direction plus itself).

Directions of Alignment

There could be more than one direction of alignment according to the alignment principle. It is possible that a cell contributes to an alignment in the 3-7 direction and in the 1-5 direction. In such cases, a non-binary CA is utilised. A cell is considered as alive if it has a value greater than 1, the state of a cell is set according to the following table:

TABLE 2 1-5 3-7 Value state 0 0 0 dead 0 1 1 alive 1 0 2 alive 1 1 3 alive

In view of this, the final cross-neighbor update rules, combining the information of Table 2 and FIGS. 10A and 10B are:

TABLE 3 c7 c5 c3 c1 c0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 1 1 0 0 1 0 0 0 0 1 0 1 1 0 1 1 0 0 0 1 1 1 2 1 0 0 0 0 1 0 0 1 0 1 0 1 0 1 1 0 1 1 2 1 1 0 0 0 1 1 0 1 1 1 1 1 0 2 1 1 1 1 3

Also the value of a cell need not be the same across all generations. For instance, a cell will have a value 3 in the first generation, but in the next generation, it will have a value 1 because its neighbors will be dead in that generation. The longest alignments in any direction are tracked using this technique. Finally, only the cells of the longest alignments will survive.

Extensions

-   1. The algorithm discussed till now considers aligning the two     structures without gaps. By a variation which involves considering     the plus-neighbors and altering the update rules to maintain a gap     cell to be live into the subsequent generations, the algorithm can     also align structures with gaps. -   2. The algorithm can be run with two sets of PPBs for one of the     protein structures. The PPB Generator 20 would skip the first Amino     Acid and generate blocks starting from the second Amino Acid so that     all possible combinations are handled. -   3. The algorithm is not only useful for detecting alignments, but     can also detect motifs dotting the CA space.

Results

An algorithm embodying the invention (presented in the Appendix, and for convenience referred to as “P-algorithm”) was tested with various sets of data. The data sets were obtained from the protein data bank of the web site of the Research Collaboratory for Structural Bioinformatics non-profit consortium.

To test the P-algorithm with similar proteins and dissimilar proteins, the DALI service for Fold Classification and Structural Alignments (FSSP) was used at the web site of the European Bioinformatics Institute, which is a non-profit academic organization that forms part of the European Molecular Biology Laboratory.

The DALI server provides a percent similarity between two protein structures based on the amino acid sequences. The similarity is a measure of the percent of similar amino acids in the final alignment. The scoring for the P-algorithm uses only the structure information and not any sequence related information.

In the P-algorithm, the best alignment is when the smaller of the two structures is completely aligned with the larger of them. The score for a P-algorithm alignment calculates the percentage of cells that are contributing to the “largest alignment”. Largest alignment here is a relative term since the P-algorithm can be configured to detect alignments of any length. The results of all the runs shown below are for alignments that are greater than 5 continuous cells. Since each cell is comprised of 5 amino acids (with an overlap of two), the length of the alignments considered for scoring are greater than 13.

If Nmax be the number of cells (Penta-Peptide Blocks) in the smaller of the two structures, and Nc be the number of cells of the final alignment, then P-algorithm Score=Nmax/Nc.

-   Run 1—Proteins “ISCU_A” having 288 C-alpha atoms and “IEUD_A” having     306 C-alpha atoms. -   Run 2—Proteins “IE6K_A” having 129 C-alpha atoms and “IM5T_A” having     124 C-alpha atoms. -   Run 3—Proteins “ID0Z_A” having 310 C-alpha atoms and “IHRK_A” having     423 C-alpha atoms. -   Run 4—Proteins “IKGS_A” having 223 C-alpha atoms and “IJFT_A” having     340 C-alpha atoms.

FIGS. 12-15 shows the visual protein structure results of running the algorithm on the different data sets, respectively for Run 1-Run 4. The following table provides a comparison.

TABLE 4 % Similarity Run No. (Dali) Nmax Nc Score 1 66 143 87 61 2 23 61 18 29 3 26 154 51 33 4 8 110 11 9

The P-algorithm compares favorably with the DALI technique.

Computer Hardware

FIG. 16 is a schematic representation of a computer system 100 of a type that is suitable for executing computer software for protein structure alignment using cellular automata. Computer software, given as psuedocode in the Appendix, executes under a suitable operating system installed on the computer system 100, and may be thought of as comprising various software code means for achieving particular steps.

The components of the computer system 100 include a computer 120, a keyboard 110 and mouse 115, and a video display 190. The computer 120 includes a processor 140, a memory 150, input/output (I/O) interfaces 160, 165, a video interface 145, and a storage device 155.

The processor 140 is a central processing unit (CPU) that executes the operating system and the computer software executing under the operating system. The memory 150 includes random access memory (RAM) and read-only memory (ROM), and is used under direction of the processor 140.

The video interface 145 is connected to video display 190 and provides video signals for display on the video display 190. User input to operate the computer 120 is provided from the keyboard 110 and mouse 115. The storage device 155 can include a disk drive or any other suitable storage medium.

Each of the components of the computer 120 is connected to an internal bus 130 that includes data, address, and control buses, to allow components of the computer 120 to communicate with each other via the bus 130.

The computer system 100 can be connected to one or more other similar computers via a input/output (I/O) interface 165 using a communication channel 185 to a network, represented as the Internet 180.

The computer software may be recorded on a portable storage medium, in which case, the computer software program is accessed by the computer system 100 from the storage device 155. Alternatively, the computer software can be accessed directly from the Internet 180 by the computer 120. In either case, a user can interact with the computer system 100 using the keyboard 110 and mouse 115 to operate the programmed computer software executing on the computer 120.

Other configurations or types of computer systems can be equally well used to execute computer software that assists in implementing the techniques described herein.

Conclusion

Various alterations and modifications can be made to the techniques and arrangements described herein, as would be apparent to one skilled in the relevant art. 

1-25. (canceled)
 26. A method for detecting alignments for protein sequences, the method comprising: generating respective data sets for two protein sequences, wherein such a data set represents protein blocks having n successive C-alpha atoms of a primary amino acid sequence; determining, for the data sets, respective matrices of distances of protein atoms from said C-alpha atoms; determining a difference matrix from said two distance matrices, wherein determining the difference matrix includes determining elements of the difference matrix, wherein such a difference matrix element is a function of a sum of absolute differences of certain elements of the respective distance matrices; and determining alignments for the protein sequences based upon the difference matrix, wherein determining the alignments includes the steps of: forming a model having a set of logical cells, the cells having respective logical states; performing iterations in which the states of the cells are selectively altered; and outputting the states to a display or a user.
 27. The method of claim 26, wherein the logical states include an alive state and a dead state, forming the model includes initially setting ones of the cells to the alive state, and performing iterations includes converting ones of the cells from the alive state to the dead state.
 28. The method of claim 27, wherein the difference matrix has a number of valued elements, each of the difference matrix elements having a respective set of neighboring ones of the elements, wherein the cells of the model form a matrix such that for each of the difference matrix elements there is a corresponding one of the cells of the model, and wherein in the initially setting of ones of the cells to the alive state such a cell is selected to be set to the alive state responsive to the value of the cell's corresponding difference matrix element relative to the values of the neighbors of the difference matrix element.
 29. The method of claim 28, wherein the cells of the model have respective sets of neighboring ones of the model cells, and wherein such a cell that was initially set to the alive state is selected in one of the iterations to remain in the alive state responsive to whether certain neighboring ones of the model cells are in the alive state.
 30. The method of claim 29, wherein the certain neighboring ones of the model cells are diagonally adjacent to the cell selected to remain in the alive state.
 31. The method of claim 26, wherein ones of said protein blocks overlap adjacent ones of the protein blocks.
 32. The method of claim 26, wherein n is a value from 1 to
 10. 33. The method of claim 26, wherein n is
 5. 